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Direct numerical simulations of 
turbulent convection with a variable 
gravity and Keplerian rotation 

By W. CABOT 


1. Motivation 

Thermal convection has been proposed as a possible mechanism for genera- 
tion and maintenance of turbulence in the inner accretion disk regime of the 
primordial solar nebula. Resulting Reynolds stresses would produce a torque on 
the disk, causing “viscous” spreading and the eventual dispersal of the gas com- 
ponent of the disk; this process was coeval with and coupled to the formation 
of the planetesimals and, ultimately, planets. Turbulent convection under solar 
nebula conditions has been described to date by ad hoc and mostly untested 
models of thermal convection and turbulence. The models of Lin & Papaloizou 
(1980) and Cabot et al (1987), in particular, give vastly different results in terms 
of the vertical distribution and intensity of turbulence, the efficiency of convec- 
tion, and structural stability properties. There are yet few, if any, conclusive 
astronomical observations of solar nebula analogues, nor does indirect evidence 
from planetary compositions provide clear constraints on these models. It is, 
therefore, of fundamental interest to design experiments with the basic physical 
features of the solar nebula in order to constrain old and new models. Although 
solar nebula conditions cannot be reproduced in the laboratory, numerical sim- 
ulations of hydrodynamic flows, which have been very successful in describing 
aerodynamic flows, can be suitably modified to provide “experimental” data for 
solar nebula modelling. 

2. Objectives 

The goals of this project are (l) to modify an extant, “proven” hydrodynamics 
code with the most important features of the solar nebula and other thin accre- 
tion disks: buoyancy terms to generate convection, internal heating representing 
the release of gravitational potential energy, a variable gravity linearly propor- 
tional to the distance from the vertical midplane due to centrifugal balance, 
rapid rotation with axis aligned with gravity, and Keplerian rotational shear; 
(2) to determine the effect that these features have on the turbulent convec- 
tion by introducing them individually and to determine the cumulative nature 
of the turbulent convection for accretion disk conditions; and (3) to model the 
convection (viz., the convective heat flux) and the turbulence (viz., turbulence 
intensities, heat dissipation, and Reynolds stresses). In this manner, prior solar 
nebula models can be tested and their deficiencies rectified. 
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3. Work to Date 

Uniform rotation. The code for direct numerical simulations of a incompress- 
ible flow in a semi-infinite channel (Kim, Moin & Moser, 1987) was modified 
to include buoyancy in the Navier-Stokes equation in the Boussinesq approxi- 
mation, the variable gravity, and uniform rotation. Both no-slip and no-stress 
boundary conditions on the vertical walls were implemented (the latter being 
more realistic for the solar nebula). Sequences of simulations with varying ro- 
tation rates were performed, from no rotation to rates approaching the limit 
of stability. A linear analysis was performed to determine theoretical critical 
rotation rates for marginal stability. It was found that the convective fluxes 
in the numerical simulations become negligible between the marginal limits of 
stationary and oscillatory convection, and there may also be finite amplitude 
effects occurring in this rotation regime. 

The turbulence intensities and efficiency of the convection were found to be 
virtually undiminished in the midchannel region where buoyancy vanishes; this 
region posed a major uncertainty in prior solar nebula modelling, and this result 
not only clarifies the nature of the turbulence there but allows some simplifi- 
cation in modelling. An extension of the standard stellar mixing length model 
for convective heat fluxes and speeds was developed to include viscosity and 
rotation and was found to agree qualitatively with the rotational dependence of 
the simulation results, although convective fluxes and speeds in the model were 
found to be a factor of 3 too low for standard values of coefficients. Canuto & 
Goldman’s (1985) turbulence model, used in Cabot et aV s (1987) solar nebula 
model, also gave a qualitatively similar dependence of convective heat flux and 
speed on rotation, but gave convective fluxes that were SO times too low. 

The vertical distribution of convective buoyancy production of turbulence ki- 
netic energy has a bimodal shape with peaks on either side of midchannel and 
vanishing at midchannel. For low rotation rates, nonlinear diffusion (primarily 
by convective transport) redistributes the turbulence so that it is dissipated with 
a flat interior profile. For high rotation rates, the overall diffusion is suppressed 
(with diffusion by pressure fluctuations now dominating) such that dissipation 
near midchannel is depressed. This result further supports the assumption made 
by Cabot et al (1987) that heat dissipation is evenly distributed throughout the 
convective region, than it does the model by Lin <fc Papaloizou (1980), in which 
the turbulence intensities and dissipation is sharply peaked at the midplane. 
Since the heat dissipation is ultimately the heat source in the solar nebula, it 
is important to assess the sensitivity of the heat dissipation distribution to that 
of the internal heat source. Tests with centrally concentrated internal heating 
show that the interior heat dissipation profile is very insensitive to that of the 
heat source, which is again consistent with the assumption that the dissipation 
will have a relatively flat interior distribution. 
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Rotational shear. Linear rotational shear is treated in the numerical integra- 
tion by working in a cosheared reference frame, in which the governing equations 
remain horizontally homogeneous, and implementing Rogallo’s (1981) remesh- 
ing transformation in the homogeneous horizontal directions, which allows one 
to circumvent the tendency of the numerical grid to become overly distorted as 
it follows the sheared flow. The remeshing transformation was implemented in 
the channel code and has been extensively tested. The tendency of the stream- 
wise length scales to become elongated and the spanwise length scales (in the 
direction of the shear) to become compressed requires balancing horizontal box 
and mesh sizes in the simulations in order to optimize the numerical resolution, 
and runs are currently under way in order to determine the optimal numerical 
domain(s). Several simulations were carried out with less than optimal but ac- 
ceptable resolution, which have provided the preliminary results subsequently 
quoted. 

Linear stability analysis predicts that the rotation profile is centrifugally sta- 
ble if the specific angular momentum gradient increases with distance from the 
rotation axis (5 > — 2H, where S is the shear rate and Cl the rotation rate). 
This limit was tested and verified with the numerical code with no buoyancy 
forces in effect. This indicates that the rotational shear flow for Keplerian ro- 
tation (5 = — 3fl/2) is centrifugally stable and, by itself, cannot maintain the 
turbulence in accretion disks. Another verification of this proposition was found 
by allowing a convecting, rotationally sheared flow simulation with well devel- 
oped Reynolds stresses due to the shear to decay when the internal heating was 
extinguished; both the convective and shear production rates decayed to zero, 
indicating that it is convection alone that drives the turbulence in this flow and 
indirectly maintains Reynolds stresses by the “catalytic” effect of the rotational 
shear. 

Preliminary results for Keplerian rotation at a moderate rotation rate with 
both no-slip and no-stress (coshearing) vertical walls show well developed 
Reynolds stresses in the horizontal components (— uuJ). The corresponding shear 
correlation coefficient was found to be 0.21 for no-stress walls and 0.23 for no- 
slip walls, compared to 0.45 for homogenous, unidirectional, plane shear flows 
(Townsend, 1956); these values were very uniform across the channel, even at 
the walls. The shear production rate of kinetic energy was about 16% of the 
convective production rate. Simulations at different rotation rates are needed to 
determine the robustness of these results. For uniform rotation with the same 
inertial frequency, the rms vertical velocity (v rms ) was found to exceed the equal 
streamwise (u rm „) and spanwise (uVm») components by about 1.4, because it is 
the component directly receiving the production by convective buoyancy. For 
the Keplerian rotational shear case, v rmt « 1.3u; rma and w rm , « 1.4u rmg , 
compared to the result for unidirectional, plane shear flow of u rmt > v rm , > 
u) rma . Because convective buoyancy is driving the turbulence in the numerical 
simulation, one expects v rms to still be the largest. The discrepancy in the 
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horizontal components may be due to rotation: The production rate in the 
streamwise component is — (2fi + S) uw — — fi/2 tluJ < 0, rather than —Suw > 0 
without rotation, and the spanwise production rate is 20 uw > 0. 


4. Future Work 

In upcoming work, numerical simulations using sequences of rotation rates 
with Keplerian shear will be constructed. The results will be used to test the 
robustness of the value of the shear correlation coefficient and the ratio of shear 
to convective production, or to suggest ways to model whatever variations arise. 
The stability behavior of the rotationally sheared flow at very rapid rotation 
rates will be compared with uniform rotation results to determine any effects of 
the shear; also the effect of shear on the modelling of convective fluxes and speeds 
will be examined. The interior distribution of kinetic energy dissipation with 
both convective and shear production will be examined for different rotation 
rates; modelling for the dissipation and the Reynolds stress will be attempted. 
Simulations with different Rayleigh and Prandtl numbers will be performed in 
order to assess their effects on the preceding key quantities with the viscosity- 
independent product of Rayleigh number and Prandtl number held close to solar 
nebula model values (~ 10 4 to 10 s ) ; lower Prandtl number simulations will allow 
higher, more realistic rotation rates to be accessed. At this stage, much more 
physical and realistic models and paxametrizations of turbulent convection in 
thin accretion disk environments, such as the solar nebula, should be attainable 
than have been available heretofore. 

In the realm of the more distant future are a number of modifications that 
would make the numerical simulations more realistic. The problem of unrealistic, 
impermeable vertical boundaries and the imposed vertical scale of the convec- 
tive region can possibly be relaxed by introducing compressibility to the code, 
at least in the form of a variable vertical density. For coverage of many density 
scale heights, the vertical boundary region would be buffered by a convective sta- 
ble region and the convective motions would select their natural vertical scale; 
however, implementation of this change will require extensive modifications of 
vertical integration scheme in the present code, probably requiring finite dif- 
ferencing rather than the present Chebyshev polynomial expansion. Since the 
dissipation of kinetic energy should properly be equated with the internal heat 
source, ways to make a self-consistent “bootstrap” will be explored; this will 
provide a self-consistent way to specify the energy source, and will indicate if 
the coupling between energy source and turbulence leads to any unanticipated 
behavior. 

This work is being carried out in collaboration with the Space Science division 
at NAS A/ Ames (J. Pollack and P.Cassen) and with NASA/GISS in New York 
(O.Hubickyj and Y.Canuto). 
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